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Abstract 

We present the first auxiliary field Monte Carlo calculations for a rare earth 
nucleus, 170 Dy. A pairing plus quadrupole Hamiltonian is used to demon- 
strate the physical properties that can be studied in this region. We calculate 
various static observables for both uncranked and cranked systems and show 
how the shape distribution evolves with temperature. We also introduce a 
discretization of the path integral that allows a more efficient Monte Carlo 
sampling. 
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I. INTRODUCTION 



Although mean-field models of rare-earth nuclei can describe the gross interplay between 
collective and single-particle degrees of freedom |T[, it is of obvious interest to pursue a 
fully microscopic description to the greatest extent possible, particularly as the angular 
momentum and/or temperature is increased. As we have previously developed and demon- 
strated Monte Carlo methods for treating very large basis shell model calculations [0|§, it 
is natural to attempt to extend them to the rare-earth region. Toward this end, we present 
here the first auxiliary field Monte Carlo calculations of observables in a rare-earth nucleus. 
For demonstration purposes, we have chosen the pairing plus quadrupole Hamiltonian and 
consider the mid-shell nucleus 170 Dy. 

As we move up to the rare earth region from the sd- and /p-shell spaces, new physical and 
computational problems must be solved. The neutrons and protons occupy different major 
shells, and so isospin symmetry is lost. Further, the numbers of single-particle orbitals and 
valence particles are larger. Most importantly, the long autocorrelations times encountered 
in Metropolis sampling [||,f|] °f continuous auxiliary fields have to be avoided by a novel 
discretization of the fields based on gaussian quadrature. The small level spacing of these 
larger nuclei also requires calculations at lower temperatures. 

Our presentation is organized as follows. In Section II we review the auxiliary field Monte 
Carlo method, discuss our choice of the Hamiltonian, and introduce a discretization of the 
path integral to facilitate the Monte Carlo sampling. We present results of static observables 
and deformation distributions for 170 Dy at various temperatures and cranking frequencies 
in Section III, and draw several conclusions about future directions. 

II. FORMALISM 

Details of the auxiliary field Monte Carlo approach to the shell model have been pre- 
sented elsewhere so that we only outline the method here. Given some many-body 
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Hamiltonian H, we desire a tractable expression for the imaginary time evolution operator, 

U = exp(-(3H) , (2.1) 

where (3 has units of inverse energy (MeV -1 , where h — 1) and (3~ x can be interpreted as 
the temperature of the system. The operator if is a generalized many-body Hamiltonian, 
and may contain terms such as — fiN in the grand-canonical ensemble, or —uJ z in cranked 
systems. (In these cases, \i is the chemical potential and u is the cranking frequency.) The 
partition function is 

Z = Trexp (-/?#) , (2.2) 

from which we can construct the thermal expectation value of an operator O as 

( 0) = I T r[Oexp (-/?#)] . (2.3) 

Here, Tr is the trace over many-body states of fixed (canonical) or all (grand-canonical) 
particle number. 

We restrict ourselves to generalized Hamiltonians that contain at most two-body terms. 
The Hamiltonian can then be written as a quadratic form in some set of convenient operators 
{O a }, 

H = J2e a O a +Wv a 2 a , (2.4) 

a * a 

where we have written the quadratic term in diagonal form. These operators are typically 
either one-particle (density) or one-quasiparticle (pairing), and the strength of the two-body 
interaction is characterized by the real numbers V a related to the two-body matrix elements. 

For H in the quadratic form ( |2.4| ), we can write the evolution operator as a path integral. 
The exponential is first split into N t 'time' slices, (3 = N t A/3, so that 

U = [exp {-A(3H)] Nt . (2.5) 

A Hubbard- Stratonovich (HS) |5|] transformation is then performed on the two-body term 
in each time slice n = 1, . . . ,N t yielding for the evolution operator 

3 



U = [exp (—A(3H)] Ni ~ J V[a]G(a)U a , (2.6) 
where the integration measure is 



2>M=nn^~(^J^Y > ( 2 - 7 ) 



n=l a 



V 2vr 



the Gaussian factor is 



G{v) = exp (- £ ~ A/3 | K | , (2.8) 

the one-body evolution operator is 

U a = U Nt U Nt _ 1 ...U l , (2.9) 

with [/„ = exp [— A/3/i CT (r„)], and the one-body Hamiltonian is 

K(r n ) =J2(e a + s a V a a an )O a , (2.10) 

where the phase factor s a is ±1 if V a < 0, and is ±i if V a > 0. Each real variable a an is the 
auxiliary field associated with O a at the time slice n. 

Expectation values of observables are calculated through 

where 

C(a) = Trf/ CT (2.12) 

is the partition function for each field configuration and 

<0>. = ™ . (2.13) 

To perform the integration via Monte Carlo methods, we introduce the non-negative 

integration weight W(a) — Q{o) \ ((a) |, and the "sign" $(cr) = C(°")/ I C( cr ) I; so that the 
expression (|2.11 ) for an observable becomes 



{U} ~ fV[a)W(vMv) ■ {2AA) 

The integral ( P-14|) can be evaluated by Monte Carlo methods using samples generated by 
the algorithm of Metropolis et al. as described in Ref. ||. However, in view of the large 
number of auxiliary fields involved (some 10 5 ), the successive field configurations are highly 
correlated for a reasonable acceptance fraction (~ 0.5), the autocorrelation length being 
over 200 sweeps. To generate uncorrelated samples more efficiently, we have approximated 
the continuous integral over each a an by a discrete sum derived from a Gaussian quadrature. 
In particular, the relation 
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A/3VO A /2 



daf{a)e AI3V ° (2.15) 



is satisfied through terms in (A/?) 2 if f(a) = | [5(cr — a ) + 5(cr + a ) + 45(cr)], where a = 
(3/VAP) 1 / 2 . (Note that commutator terms render the HS transformation accurate only 
through order A/3 anyway.) 

In this way, each a an becomes a 3-state variable and the path integral becomes a (very 
large) path sum, which can be sampled using the algorithm of Metropolis et al. We have 
found that this discretization reduces the correlation length to only five sweeps, and thus 
increases our efficiency by a factor of 40 relative to the continuous case, with no loss of 
accuracy. 

To describe rare-earth nuclei, we choose the Z = 50-82 shell for protons (2sl<i0g7/20/iii/2) 
and the iV = 82-126 shell for neutrons (2pl /Ogg/2 0^13/2) . The Hamiltonian we have chosen 
is of the pairing plus quadrupole form given by 

H = E e a a\fl a - g v P\P v - g n P\P n - f Q p • Q n , (2.16) 

where e a are the single particle energies, a£, and a a are the anti-commuting creation 
and annihilation operators associated with the single particle state a, Pp( n )->Pp(n) are the 
monopole pair creation and annihilation operators for protons and neutrons, and Q p ( n ) is 
the quadrupole-moment operator. The single particle states a are defined by the complete 



set of quantum numbers nljmt z , denoting the principal, orbital angular momentum, total 
single-particle angular momentum, ^-projection of j, and the third component of the isospin 
quantum numbers, respectively. Single-particle operators are thus represented by matrices of 
dimension 32 and 44, for protons and neutrons, respectively. The pairing strengths g p ( n ), the 
quadrupole interaction strength x, and the single particle energies used in our calculations 
are given in Table | (taken from Ref. ||). 



III. THE CALCULATION AND RESULTS 

To demonstrate the power of our methods, we study the mid-shell nucleus 170 Dy (16 
valence protons and 22 valence neutrons), which requires some 10 21 m-scheme determinants. 
We have used Aft = 0.0625 and N t = 8 to 64 time slices. 

In addition to grand-canonical calculations, we performed canonical analyses of the fields 
generated through grand-canonical sampling. In particular, canonical observables (subscript 
"c" ) are given in terms of the grand-canonical sampling (subscript "<?" ) as 

(n) SV[cr]w g (*)* e (<T) [U°)IU°)W)°c , n 

[ )c JV[a]W 9 (a)<S> c (a)[( c (a)/( g (a)} ' 1 ' > 

The fluctuations in [( c (&) / Cg( (J )} determine how precise this evaluation can be. We find 
that the fluctuations are less than 10% (as the number distribution in the grand-canonical 
ensemble is small), so that canonical observables can be calculated with good precision. 

In Fig. ([I]) we show the static observables for the uncranked system in both the grand- 
canonical and canonical formalisms. We calculated observables canonically, using grand 
canonical fields, up to (3 = 2.0 in order to demonstrate that for these nuclear systems, either 
method may be used in this kind of calculation. Note that as the temperature decreases, 
the nucleus becomes more deformed. Relaxation of the expectation values of H and J 2 is 
also clearly seen. The sign in these cases is identically one for all temperatures. 

Cranking calculations in which H — > H — ujJ z have also been performed. The systematics 
are shown in Fig. |2], where we display ($), (H), (Q 2 ), (J 2 ), —g(P^P), and the moment of 
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inertia, (X 2 ). Note that the sign degrades quite rapidly with increasing u, making cranking 
calculations at lower temperatures difficult. Moments of inertia were calculated from X 2 = 
d(J z )/du = (3[{Jl) — (J z ) 2 ]- At high temperatures, the nucleus is unpaired and the moment 
of inertia decreases as the system is cranked. However, for lower temperatures when the 
nucleons are paired, the moment of inertia initially increases as we begin to crank, but then 
decreases at larger cranking frequencies as pairs break; Figure 2 shows that the pairing gap 
also decreases as a function of u. It is well known that the moment of inertia depends on 
the pairing gap 0, and that initially X2 should increase with increasing u. Once the pairs 
have been broken, the moment of inertia decreases. These features are evident in the figure. 

In addition to static observables, we have calculated the nuclear shape distributions. 
These distributions give a clear description of various shape and phase transition phenomena 
||. To obtain a detailed picture of the deformation, we use the components of the quadrupole 
operator = r 2 !^* . Rotational invariance of the Hamiltonian implies that the expectation 
value of each component of vanishes on average; however, each Monte Carlo sample will 
have finite values, from which we construct the quadrupole tensor = 3xiXj — S^r 2 . 
The eigenvalues of the latter tensor then lead directly to the deformation parameters ||. 

Figure |3| shows the evolution of the shape distribution for 170 Dy at inverse temperatures 
T _1 = 0.5, 1.0, 2.0, and 3.0 MeV -1 . These contour plots show the free energy 7), 
obtained from the shape probability distribution, P(f3, r y), by 

F(/?, 7 ) = -Tln^l^, 
p 6 sin 37 

where the /3 3 sin 37 is the metric in the usual deformation coordinates. As is seen from the 
plots, deformation clearly sets in with decreasing temperature. At high temperatures, the 
system is nearly spherical, whereas at lower temperatures, especially at T _1 = 3.0, there is 
a prolate minimum on the 7 = axis. 

These calculations demonstrate how auxiliary field Monte Carlo methods can be extended 
to rare earth nuclei. For the first time we have used different major shells for protons and 
neutrons in such a calculation. We have demonstrated how to obtain canonical information 



from grand-canonical sampling, and have introduced a discretization of the field integrals 
that allows for a much more efficient sampling of the integrand. Results for 170 Dy with a 
pairing + quadrupole Hamiltonian show qualitatively the expected behavior as a function 
of both cranking frequency and temperature. While the schematic nature of the Hamilto- 
nian precludes reading too much into our results, such calculations could be used to check 
various approximation schemes. Furthermore, the recently demonstrated ability to handle 
more realistic shell-model hamiltonians via Monte Carlo methods flu]] should allow these 
calculations to move beyond the schematic level. 
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FIGURES 

FIG. 1. Canonical and grand-canonical observables for the uncranked 170 Dy system as a 
function of (3. Circles represent results of canonical projection of the grand-canonical fields, 
while squares show the grand-canonical results. Lines are drawn to guide the eye and all er- 
ror bars are smaller than the symbol sizes. We show expectation values of the energy (H), 
the isoscalar quadrupole moment (Q 2 ), the valence nucleon number (N) and number vari- 
ance (AN) = \J (N 2 ) — (N) 2 (grand-canonical only), the squared angular momentum (J 2 ), and 
—g(P^P), the expectation value of the pairing terms in the hamiltonian. 

FIG. 2. Grand-canonical observables for 170 Dy at various cranking frequencies and tempera- 
tures. We show the average sign (<&), the isoscalar quadrupole moment (Q 2 ), the energy (H), the 
square of the angular momentum (J 2 ), the moment of inertia (I2), and the expectation value of 
the pairing terms in the hamiltonian —g(P^P). Error bars where not shown are approximately the 
size of the symbols, and lines are drawn to guide the eye. 

FIG. 3. Contours of the free energy (as described in the text) in the polar-coordinate plane 
for the 170 Dy system. Inverse temperatures are from 0.5 (top) 1.0, 2.0, and 3.0 MeV -1 (bottom). 
Contours are shown at 0.3 MeV intervals, with lighter shades indicating the more probable nuclear 
shapes. As our calculations become indeterminate at 7 = 0, we have truncated these plots at small 
7- 
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TABLES 



TABLE I. Physical parameters used in these calculations. 





Protons 




Neutrons 


Orbital 


e (MeV) 


Orbital 


e (MeV) 




-2.24 


0^9/2 


-3.540 


14/2 


-1.979 


L/7/2 


-3.211 


0/in/2 


-1.120 


0«l3/2 


-2.240 


14/2 


-0.122 


2p 3 /2 


-1.20 


2si/ 2 


0.000 


1/5/2 


-0.933 






2P1/2 


0.00 


5p = 


+0.160 MeV 


5n 


= +0.131 MeV 



X = +0.054 MeV fm 4 
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